#ifndef SAVE_H
#define SAVE_H

#include <functional>
#include <deque>
#include <assert.h>

class SavedFunction : 
    public unary_function<double, double>
{
  private:
    double (*f)(double);
    double step;
    std::map<double,double> fmap;
    typedef std::map<double,double>::iterator mapit;

  public:

    SavedFunction(double (*_f)(double), double _step) : f(_f), step(_step) {}

    double operator()(double x)
    {
      if (fmap.size() == 0) {
	double fx = (*f)(x);
	fmap[x] = fx;
	return fx;
      } else {
	mapit it1 = fmap.upper_bound(x);
	while (it1 == fmap.begin()) {
	  double newx = it1->first - step;
	  double newfx = (*f)(newx);
	  fmap[newx] = newfx;
	  it1 = fmap.begin();
	  if (it1->first <= x) it1++;
	}
	mapit it0 = it1; it0--;
	while (it1 == fmap.end()) {
	  double newx = it0->first + step;
	  double newfx = (*f)(newx);
	  fmap[newx] = newfx;
	  it0 = it1 = fmap.end(); it0--;
	  if (it0->first > x) {it0--; it1--;}
	}

	assert(it0->first <= x);
	assert(it1->first > x);

	double x0 = it0->first;
	double x1 = it1->first;
	double f0 = it0->second;
	double f1 = it1->second;
	double fx = (x-x0)/(x1-x0)*(f1-f0) + f0;
	return fx;
      }
    }
};

#endif
